Journal of Mathematical Biology

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [&ge;] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

We show that the p53-Mdm2 gene interaction can be understood as a form of antithetic feedback, but with mixed feedback loops. In particular, we show that this viewpoint sheds light on the emergence and robustness properties of p53's oscillatory behaviour, as observed experimentally.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Species at the highest trophic levels are generally the least resilient to perturbations, yet the factors that determine their persistence remain poorly understood. We ask how the intrinsic self-regulation of a basal host governs the persistence of the trophic levels above it. To address this question, we use a minimal discrete-time model of a host-parasitoid-hyperparasitoid system motivated by the biology of trigonalid wasps, obligate hyperparasitoids that complete development only when their caterpillar host is also attacked by a primary parasitoid. We compare two forms of host self-regulation in the absence of natural enemies: scramble (Ricker) fecundity, appropriate when the parasitoid is the host's principal enemy and which recovers the classical model of Beddington and Hammond as a special case; and contest (Beverton-Holt) fecundity, appropriate when the host is held in check by many additional enemies. Using numerical continuation, Lyapunov exponents, and direct computation of the model's trajectories, we find that the form of host regulation is decisive for persistence at the top trophic level. Under scramble host dynamics, coexistence is extremely fragile with respect to both parameter values and initial conditions. Under contest host dynamics, coexistence is robust across broad parameter ranges and loses stability only gradually, through a supercritical Neimark--Sacker bifurcation and a subsequent torus-doubling route to chaos. Therefore, in at least some trophic networks, the resilience of the highest levels can be determined at the base of the network.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Ecological coarse-graining of individuals or species into aggregate variables, such as total biomass or regional abundance, is common practice, yet differences among species can still influence aggregate dynamics. We propose a metric for quantifying this irreducible dynamical heterogeneity. For small fluctuations near a stable equilibrium, we compare the aggregate fluctuation spectrum of the full dynamical system with that of the best Ornstein-Uhlenbeck model having the same dimension as the aggregate. The remaining spectral mismatch quantifies the dynamical effects of hidden differences that the fitted aggregate model cannot capture. For weakly interacting random community models at fixed collective interaction strength, we show that increasing species richness can either reduce or increase this heterogeneity, depending on how species interact. Our framework helps understand when ecological and evolutionary differences average out versus when their dynamical effects persist.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Devireddy et al. [1] previously showed that an enthalpy-based, macroscale model simulated tissue freezing histories are in close agreement with a coupled model incorporating cellular water transport and intracellular ice formation (IIF). That comparison, however, was performed for a limited set of tissue geometries, cooling conditions, and biophysical parameters. Here, we revisit that conclusion by reimplementing the published model from its governing equations and systematically expanding the parameter space to determine when microscale processes alter the predicted macroscale freezing response. We incorporate both surface catalyzed nucleation (SCN) and volume catalyzed nucleation (VCN) of intracellular ice with cellular osmotic dehydration and examine a broad range of tissue dimensions, cooling rates, convective boundary conditions, membrane permeabilities, cell sizes, activation energies, intracellular water fractions, and nucleation rates. Linearization of the cellular water transport equation about osmotic equilibrium yields an osmotic relaxation timethat can be compared with the local residence time of tissue undergoing phase change. Their ratio defines a dimensionless group, that organizes the transition between regimes in which the uncoupled enthalpy formulation and the coupled micro and macroscale formulation produce similar or substantially different results. Across the one-dimensional parameter space examined here, the two models remain within a prescribed thermal-error bound for the dimensionless group below approximately 0.1, whereas the discrepancy increases systematically above this threshold. Model divergence also depends on intracellular ice nucleation, i.e., when nucleation is sufficiently slow, cellular water transport can become rate limiting for latent heat release. Conversely, sufficiently rapid intracellular ice formation reduces the macroscopic discrepancy by providing an additional pathway for latent-heat release. Multidimensional simulations further show that geometry alone can move a nominally conventional freezing protocol into the coupled regime. In particular, multidirectional cooling can reduce the local residence time sufficiently to produce substantial differences in predicted frozen volume even when the corresponding thermal histories appear similar. Finally, we formulate an inverse approach for estimating IIF parameters from tissue scale thermal histories. Synthetic simulations indicate that parameter recovery is feasible within a finite range of cooling conditions, providing a testable strategy for future experimental measurements. Together, these results provide a mechanistic criterion for deciding when microscale biophysical processes must be retained in macroscale tissue-freezing models.

Game theory in biology started as an attempt to model animal contests by making assumptions about the fitness costs individuals pay in aggressive interactions. Among the different approaches, a notable one is to construct mechanistic models, using assumptions about behaviours and cognitive processes. Reinforcement learning is an important cognitive process, and we use it here to model pairwise and coalitional contests. We study situations where territorial neighbours become so-called dear enemies. A defending individual can get help from a neighbour through a defender-neighbour coalition against a challenger attempting territory takeover. In our model, coalition members have an advantage in contests in terms of costs of aggression, with the challenger being exposed to aggression from both members and each member only receiving part of the aggression from the challenger. We find that notable costs of border conflicts between neighbours, together with substantial advantages for coalitions, favour intervention. Neighbours intervene when they can help weaker defenders, in this way avoiding costs of border renegotiation with a new and potentially stronger territory owner. We also introduce different forms of perceived costs into the model, which we show correspond to the range from pure mutual assessment to partial self-assessment.

Haldane's principle relates equilibrium mutation load to deleterious mutational input rather than to selection strength. We study this balance in a homeostatically regulated, damage-structured population where mutation occurs during reproduction, repair restores individuals toward lower-damage states, and nonlinear demographic feedback maintains constant population size. An exact demographic identity shows that load at every time equals excess reproductive effort weighted by mean relative fecundity. Without repair, the equilibrium branch with positive pristine-class frequency is exactly solvable, and along it the load is exactly the mutation probability times the realized per-capita reproductive turnover rate; when the terminal cap is removed, this is precisely the realized deleterious input, for every mutation probability below one. For rare mutation it reduces to the mutation probability times the pristine turnover rate, with the local fitness contrast canceling at leading order. With a finite terminal class the branch instead reaches a closure-dependent boundary collision. Repair modifies the rare-mutation balance by partitioning first-order removal from the first damaged class between selection and restoration, multiplying the leading load by the selective share of that removal. An exact binary solution extends these results to finite mutation pressure: the unrepaired interior equilibrium collides nonhyperbolically with the fully damaged boundary and converges algebraically at the threshold, whereas any positive repair rate restores a unique interior equilibrium for positive mutation; at the no-repair threshold, its displacement from the boundary scales as the square root of the repair rate as repair vanishes. These results separate the familiar competing-transition structure of mutation--selection balance from the specifically homeostatic effects of turnover-dependent mutation input and nonlinear demographic feedback.

We develop a general geometrical framework for growth processes in which the growth rate depends exponentially on the number of active growth-supporting components. The framework separates the fraction of active components from the size of the accessible component space, allowing qualitatively different growth laws to arise from the same underlying construction. Bacterial population growth provides an example of a fixed-size component space controlled primarily by activation and deactivation, whereas examples such as human population growth or speciation provide an expanding component space in which interactions between individuals promote further growth. This motivates a Replicative geometrical law, dX/dt= X exp({beta} X^{gamma}), whose exponent characterizes the scaling of the effective growth-supporting component space and is tested against historical human-population regimes and eukaryotic lineage-through-time data. For innovation indicators such as patents and arXiv manuscript submissions, where an existing item need not itself reproduce, the same construction motivates a Direct Geometrical production law without the prefactor X. The framework therefore provides a common description of fixed and expanding component spaces while distinguishing Replicative from Direct Geometrical growth.

Kinetic parameters in engineered biomolecular circuits are typically known only within intervals and are not precisely tunable; controller synthesis is therefore naturally posed over sets of gains rather than at a single nominal design. We formulate robust controller tuning under parametric plant uncertainty as a universally quantified set-inversion problem: given a plant-parameter box, compute the controller gains that stabilize every plant in the box. Robust stability is decided by Hurwitz sign conditions, which are bounded over the joint gain-parameter domain using validated range enclosures, yielding machine-checkable certificates of local exponential stability of the deterministic equilibrium. Numerical experiments show that the certified region depends strongly on the range enclosure. For an antithetic integral feedback (AIF) controller acting on a two-species gene-expression cascade, a tensor-product Bernstein enclosure certifies 83.2% of the design plane at bisection tolerance epsilon=0.05 with 324 terminal boxes, whereas four baseline enclosures certify between 50.9% and 58.9% on the same paving. For a fifth-order loop, the certified fraction increases from 34.4% to 73.8%, and the certified and excluded terminal boxes confine the robustly stabilizing area to a window of about 2.1 percentage points. Finally, for a refining variant of the algorithm in which the plant box is refined alongside the design box, the undecided area satisfies |U(epsilon)|=O(epsilon) under a non-degeneracy assumption on the stability boundary.

Allergic conjunctivitis is an inflammation of the conjunctiva caused by allergen; it is common disorder affecting up to 40% of the population. In this work, we study seasonal allergic conjunctivitis (SAC), also called "hay fever eyes", which is caused by exposure to airborne pollens. We develop a mathematical model quantifying the ocular immune system response to the allergens. First, we present a simplified qualitative description of the immunopathogenesis of SAC. Then, we express each chosen immunopathological mechanism mathematically to construct a system of thirty-one ordinary differential equations. We compare summary statistics of the predicted observable immune signals to experimental measurements and find our model captures key qualitative features of SAC progression. We then compare our predicted time series of histamine concentration to symptom scores and find a strong correlation suggesting the model predicts relevant clinically trends. Next, we calibrate the model through multi-step process. We find the most influential parameters are the production and depletion rates of IL-4, and the production rates of IL-5 and IL-8. These cytokines are targeted in treatments for asthma, atopic dermatitis, and severe eosinophilic associated disorder and suggest potential therapeutic targets for SAC. Our calibrated model mimics most of the summary statistics of the experimentally observable immune signals with discrepancies for IL-5 and IL-13 indicating that additional immunopathological mechanisms could be important.

Oncogene amplification on circular extrachromosomal DNA (ecDNA) has been linked to poor prognosis and higher treatment resistance in multiple types of human cancer. ecDNA are mobile genetic elements lacking centromeres that are partitioned unevenly into daughter cells at mitosis. While random segregation of ecDNA contributes to gene-copy-number heterogeneity among tumour cells, how ecDNA contribute to phenotypic heterogeneity, via mRNA and proteins, and how this affects targeted therapy outcomes is still not understood. In fact, cancer cell populations with ecDNA show remarkable heterogeneity in mRNA concentrations, which is not explained by gene-copy-number heterogeneity alone. In this paper, we use a stochastic model of gene transcription to mathematically model mRNA production and decay in cancer cells with ecDNA, extending a method that has been used to study phenotypic heterogeneity in eurokaryotic and prokaryotic cells. Using our model, we explore the space of possible mRNA distributions from a given ecDNA distribution, under different model parameter regimes. Our model is able to produce realistic mRNA distributions and we use it to estimate transcriptional parameters from single-cell measurements of ecDNA and mRNA counts in cell-line data. We discover a remarkable similarity in estimated parameters across cell lines and amplification types, indicating similar processes underlying transcription in ecDNA and chromosomal amplifications. We also model genotoxic therapy and use our model to predict the mRNA composition of persisting cells, showing the usefulness of such a modelling framework to understanding the origins of treatment resistance in ecDNA-containing cancers.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Regular environmental cycles such as day and night or the lunar cycle are ubiquitous and drive many biological processes. Many organisms have evolved endogenous time-keeping mechanisms to anticipate these cycles to accurately time their behavior or synchronize physiological processes. The evolution of such synchronization to environmental cycles and their variability is not very well understood. As an example, we use the marine insect Clunio (Diptera, Chironomidae), which synchronize adult emergence and reproduction with the lunar cycle, and show different emergence patterns across populations. We use a model constructed with replicator-mutator dynamics to provide a novel outlook at the evolution of these emergence patterns. In particular, we investigate the effect of two traits of individual phenotypic variability of emergence, the spread and shape, on population-wide emergence patterns. We find that the semilunar and lunar emergence patterns are alternatively stable and occur when the phenotypic distribution has a narrow and is roughly flat. Arrhythmic emergence occurs when the phenotypic distribution is wide and nearly flat. The semi-lunar pattern can only be observed when the parameter of rhythmicity is changed. With our model, we introduce the concepts of frequency-dependent selection and individual phenotypic variability into the theory of chronobiology, which gives us insights into the evolution of rhythmicity in these organisms.

Understanding how species interactions, ecological drift, and immigration jointly shape community composition is a central goal of community ecology. These processes determine not only average community composition, but also how often different community states occur. Here, we ask how competition reshapes the neutral distribution of community composition generated by drift and immigration. We analyze a stochastic two-species competition model with immigration and derive the stationary distribution of relative species composition near neutrality. The stationary distribution shows when communities are usually balanced and when they are instead dominated by one species. Species dominance can become common even when the corresponding deterministic model predicts stable coexistence, particularly in smaller communities. Balanced and dominance states can also occur together, producing a trimodal distribution. Competitive asymmetry further biases the distribution toward the competitively favored species. Thus, the interplay between competition, ecological drift, and immigration determines not only how much community composition varies, but also which community states are most likely to be observed. Our results provide an analytical bridge between neutral and competition theory and a framework for interpreting variation in community composition.

A major transition in the evolution of eukaryotic cells was phagocytosis: the ability of one cell to engulf another. In some cases, the engulfed cell persists, forming an endosymbiosis which is often mutually beneficial. However, it is unclear whether the interaction between the ancestors was also cooperative. There are clear benefits to developing the ability to eat a competitor, but can it ever be beneficial to eat a cooperator? We use a dynamical systems framework to identify the ecological conditions under which phagocytosis can evolve from a system of cooperating microbes.Beginning with a system of two cooperating species, we use invasion analysis to identify the parameter regimes under which a phagocytosing mutant can invade and the subsequent population dynamics. We identify a critical threshold that determines whether an invading predator will coexist with the cooperative community or replace members. Repeatedly invading with different predators, we find that the system saturates when there are three coexisting species, though this community does not have a single optimal composition. However, if we incorporate additional costs, we find there is a single evolutionarily stable community of two cooperators and a pure predator. Ultimately, we show that it is possible for phagocytosis to evolve from a cooperative community and, moreover, coexist with the original cooperating species.